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Abstract. In this paper we provide a reconstruction algorithm for piecewise-smooth functions from their 
Fourier coefficients, posessing the maximal possible asymptotic rate of convergence - including the positions 
of the discontinuities and the pointwise values of the function. This gives a positive answer to the so-called 
Eckhoff's conjecture. The vital difference of our algorithm compared to the original Eckhoff's method is 
that instead of consecutive Fourier coefficients, one should take ones that are evenly spaced throughout the 
whole sampling range. 

1. Introduction 

Consider the problem of reconstructing a function / : [— 7r,7r] — > R from a finite number of its Fourier 
coefficients 

c k {f) = ^ f(t)e- M dt, k = 0,l,...M. 
It is well-known that for periodic smooth functions, the truncated Fourier series 

M 

Sm(/)= E c h(f)e lkx 

|fc|=0 

converges to / very fast, subsequently making Fourier analysis very attractive in a vast number of applica- 
tions. We have by the classical Lebesgue lemma (see e.g. [31]) that 

max |/(a:) - £ M (/) Or) I < (3 + In M) ■ E M (/) 

— 7Z<X<7T 

where Em (/) is the error of the best uniform approximation to / by trigonometric polynomials of degree at 
most M. This quantity, in turn, depends on the smoothness of the function. In particular: 

(1) If / is (i-times continuously differentiable (including at the endpoints) and — then (see 
[401 Vol.1, Chapter 3, Theorem 13.6]) 

(1.1) E M (/) < C d ■ R ■ M- d . 

(2) If / is analytic, then by classical results of S.Bernstein (see e.g. [22 Chapter IX]) there exist constants 
C and q < 1 such that 

(1.2) E M {f)<C-q M . 

Yet many realistic phenomena exhibit discontinuities, in which case the unknown function / is only piecewise- 
smooth. As a result, the trigonometric polynomial 3\M (/) no longer provides a good approximation to / 
due to the slow convergence of the Fourier series (one of the manifestations of this fact is commonly known 
as the "Gibbs phenomenon"). It has very serious implications, for example when using spectral methods 
to calculate solutions of PDEs with shocks. Therefore an important question arises: "Can such piecewise- 
smooth functions be reconstructed from their Fourier measurements, with accuracy which is comparable to 
the 'classical' one (such as (|1.1[) or (|1.2p )"? 
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This problem has received much attention, especially in the last few decades ([21 131 IH EH [HI H21 EH E3 
[T3[T![T![2ll[2![M[ia[M[23 would be only a partial list). It has long been 

known that the key problem for Fourier series acceleration is the detection of the shock locations. Applying 
elementary considerations we have the following fact (see prrof in Appendix [X| . 

Proposition 1. Let f be piecewise d-smooth. Then no deterministic algorithm can restore the locations of 
the discontinuities from the first M Fourier coefficients with accuracy which is asymptotically higher than 
M- d ~ 2 . 

Until now, the question of whether this maximal accuracy is achievable in practice remained open. In 
our previous work |5j we have provided an explicit reconstruction algorithm (see below) which restored the 
jump locations (and subsequently the pointwise values of the function between the jumps) with "half the 
maximal accuracy. In the present paper we modify the method of [9] (see Algorithm[3| so that full asymptotic 
accuracy is achieved (Theorem Thus we settle this almost 20-year old "Eckhoff's conjecture". The vital 
difference of the new algorithm compared to the original Eckhoff's method (and its modification from [5]) is 
that instead of consecutive Fourier coefficients, one should take ones that are evenly spaced throughout the 
whole sampling range (thus we call the new method "decimated Eckhoff's algorithm"). 

We describe the general approach, as well as our previous results obtained in [9], in Section [2j The 
modified algorithm is provided in Section^ and its accuracy is analyzed in Section^ Finally, some numerical 
experiments are presented in Section [5] 

We would like to thank B.Adcock and Y.Yomdin for useful discussions. 

f .1. Notation. The space of piecewise C d functions on [— tt, it] with K jump points is denoted by PC (d, K). 

2. Previous work 

2.1. Eckhoff's method. Let us first briefly describe what has become known as the Eckhoff's method for 
nonlinear Fourier reconstruction of piecewise-smooth functions [161 1171 ITS] . 

Assume that / has K > jump discontinuities {£j},- =1 (they can be located also at ±7r, but not necessarily 
so). Furthermore, we assume that / € C d in every segment (£,-_x,£,-), and we denote the associated jump 
magnitudes at £j by 

au 3 =f {l \it)-f {l) ttj)- 

We write the piecewise smooth / as the sum / = >!' + $, where ^>{x) is smooth and periodic and $(x) 
is a piecewise polynomial of degree d, uniquely determined by , {a<i,j} such that it "absorbs" all the 
discontinuities of / and its first d derivatives. This idea is very old and goes back at least to A.N.Krylov 
([21 [27]). Eckhoff derives the following explicit representation for Q(x): 

K d 

<^) = EE a <.^o) 

(2.1) 3 = 1 1=0 

Vn {x; & = - (^TTji Bn+1 [~^r) 6 - x - ^ + 2n 

where V n (x;£j) is understood to be periodically extended to [— 7r,7r] and B n [x) is the n-th Bernoulli poly- 
nomial. Elementary integration by parts gives the following formula. 

Proposition 2. Let be given by (|2.ip . For definiteness, let us assume that co($) = f^<f>(x)dx = 0. 

Then 

K d 

(2.2) c fe (<I>) = 5>fc)-'-V,- 

n j=i i=o 

Eckhoff observed that if ^ is sufficiently smooth, then the contribution of Cfc( x P) to Ck(f) is negligible for 
large k, and therefore one can hope to reconstruct the unknown parameters {£j) a z,j} from the perturbed 
equations (|2.2|) . where the left-hand side reads Ck (/) ~ Ck ($) and k ^> 1. His proposed method was to 
construct from the known values 



{ Ck (/)} k = M - (d + 1) K + 1, M - (d + 1) K + 2, . . . , M 
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an algebraic equation satisfied by the jump points {£1, . . . , £k}, and solve this equation numerically. Based on 
some explicit computations for d = 1, 2; K = 1 and large number of numerical experiments, he conjectured 
that his method would reconstruct the jump locations with accuracy M^ d ^ 1 . 



2.2. Half-order reconstruction. In [9\ we proposed a reconstruction method based on the original Eck- 
hoff's procedure as follows. 



Algorithm 1 Half-order algorithm, [9]. 

Let / € PC (d, A), and assume that / = <£>( d ) + ^> where & d > is the piecewise polynomial absorbing all 
discontinuities of /, and iteCf Assume in addition the following a-priori bounds: 

• Minimal separation distance between the jumps 

mp|&-&| > J > 0; 

• Upper bound on jump magnitudes 

\ai,j\ < A < oo] 

• Lower bound on the value of the lowest-order jump 

\ao,j\ >B>0; 

• Upper bound on the size of the Fourier coefficients of ^: 

\c k m\<R-k- d -\ 

Let us be given the first 3M Fourier coefficients of / for M > M (d, A, J, A, B, R) (a quantity which is 
computable). The reconstruction is as follows. 

(1) Obtain first-order approximations to the jump locations {£i, . . . , £k} by Prony's method (Eckhoff's 
method of order 0). 

(2) Localize each discontinuity £j by calculating the first M Fourier coefficients of the function fj = f -hj 
where hj is a C°° bump function satisfying 

(a) hj = on the complement of — J, £j + J] ; 

(b) ftiElon + 

(3) Fix the reconstruction order d\ < [|J. For each j — 1,2,..., A, recover the parameters 

do j, . . . , ad.!,]} from the d\ + 2 equations 

c k (fj) = y e~^ k ~~h^i k = M - dt - 1,M - di, . . . ,M 

by Eckhoff's method for one jump (in this case we get a single polynomial equation (£/) = o| 
of degree di). 

(4) From the previous steps we obtained approximate values for the parameters ■fCj'j an d {ai.j}. The 
final approximation is taken to be 

e f*(/)-^E^*i;7|fe}e*+f;i;2M^^). 

ifei<M I j=i ;=o v ; j=i ;=o 



We have also shown that this method achieves the following accuracy. 
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Theorem 3 ([5]). Let f 6 PC(d,K) and let f be the approximation of order d\ < [41 computed by 
Algorithm^ T/ierfl 



< d (d, di.if, J,A,B,R) ■ M 



-di-2. 



(2.3) 



a/ j I < C 2 (d, di,if, J, A, B, R) ■ M^ 1 " 1 , I = 0, 1, . . . , d x ; 



f(x)-f (x) < C 3 (d, d lt K, J, A, B, R) ■M- d '-\ 



The non-trivial part of the proof of this result was to analyze in detail the polynomial equation p (£j ) = 
in step 3 of Algorithm [T] It turned out that additional orders of smoothness (namely, between d\ and d) 
produced an error term which, when substituted into the polynomial p, resulted in unexpected cancellations 
due to which the root £j was perturbed only by O (M~ dl ). This phenomenon was first noticed by EckhofF 
himself in [17] for d = 1, but at the time its full significance was not realized (see Section [6]). 



3. The decimated Eckhoff algorithm 

We begin the exposition by pointing out an important property of Algorithm [TJ namely that its final 
asymptotic convergence order essentially depends on the accuracy of step 3. It is sufficient therefore to 
replace this step with another method which achieves full accuracy (i.e. ~ M~ d ~ 2 ) in order to obtain the 
overall reconstruction with full accuracy. It turns out that taking instead of consecutive Fourier samples 

k = M-d-l,M-d,...,M 

the "decimated" section 



k = 



M 



(d + 2) 



• 2 



M 



(d + 2) 



..,M 



does the job. 

In what follows, we assume that the jumps locations £j are already known with accuracy at least O [M _2 J . 
Instead of the index M we shall be using 

M 



N 



del' 



.(d + 2). 

throughout and in the end we'll just switch back to M . 

Denote the single jump point by £ € [— ir, tt], and let to = e~ 1 ^ . 

The purpose is to recover the jump point uj and the jump magnitudes {do, . . . , ad} from the noisy mea- 
surements 

,k 



Ck 



( zfc ) 



J+l 



+e k , k = N,2N,...,(d+2)N, \e k \<R-t 



-d-2 



We multiply by (2tt) (ik) d+l and let otj = i a+1 ■'aa—j to get 



- ,d+l-j. 



(3.1) m k d = 2tt {ik) d+1 c k = Lu k 



ctjk? +4, 



k = N, 2N, . . . , (d + 2) N, \5 k \<R-k" 



Let z = u N . We have 



Definition 4. Let 



def 

= m k 



777 I 



d i \ def 



d+1 



i=0 



d+1 



3=0 V J 



d+l-j 



The last (pointwise) bound holds on "jump-free" regions. 
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Proposition 5. The point u = z is a root of p% (u). 
Proof. Combinatorial exercise. 

Definition 6. Let V$ denote the (d + 1) x (d + 1) matrix 

"1 N N 2 



xrd 4£ f 

v N — 



1 



2N 



(2Nf 



I (d+l)N {{d+l)NY . 
Proposition 7. The vector of exact magnitudes {ay} satisfies 

-N 



N d 
(2N)° 



((d+l)N)° 



(3.2) 



-2N 



m (d+ i)jvw 



-{d+l)N 





a 




Cil 






pt d _ 



□ 



Algorithm 2 Recovery of single jump parameters 

Let there be given the first iV 3> 1 Fourier coefficients of the function fj, and assume that the jump 
position £ is already known with accuracy O (7V~ 2 ). 

(1) Construct the polynomial 

d+l 

from the given noisy measurements mjv, m,2N, ■ ■ ■ , fh(d+2)N (|3.1|) . 

(2) Find the root z which is closest to the unit circle (in fact any root will do, see Remark [TBI below) . 

(3) Take uj = \fz- Note that in general there are N possible values on the unit circle, but since we 
already know the approximate location of u) the correct value can be chosen consistently. 

(4) Recover £ = — argw. 

(5) To recover the magnitudes, solve the linear system (|3.2[) . 



( d+ ] Jrn (j+1)N u d+1 K 



To see that there are N possible solutions, notice that 

£N-t = 2im 

t 2?r 

£ = 1 n, 

s N N 



Algorithm 3 Full accuracy Fourier approximation 

Let / £ PC (d, K) , and assume that / = $W + * where is the piecewise polynomial absorbing all 
discontinuities of /, and "J/ € C d . Assume the a-priori bounds as in Algorithm [TJ 

(1) Using Algorithm [TJ obtain approximate values of the jumps (up to accuracy ) and the 

Fourier coefficients of the functions fj, 

(2) Use Algorithm [2] to further improve the accuracy of reconstruction. 
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4. Accuracy analysis 

Theorem 8. Algorithm^ recovers the parameters of a single jump from the given noisy measurements (|3.ip 
with the following accuracy: 

R 



< Ci-N- d - 



ISj-ajl < CsRfl + ^jN-j- 1 , j = 0,l,...,d. 

An immediate consequence is the resolution of Eckhoff 's conjecture. 
Theorem 9. Let f G PC (d,K) and let f be the approximation of order d computed by Algorithm^ Then 



< C 6 (d, K, J, A, B, R) ■ M 



-d-2. 



(4.1) \ay-aijl <C 7 (d,K,J,A,B,R)-M l - d -\ I = 0,1,.. .,d; 

f(x) -f(x)\<C s (d, K, J, A, B, R) • M^" 1 . 

The remainder of this section is devoted to proving Theorem [5J 
Let us first define an auxiliary polynomial sequence. 

Definition 10. For all i, d nonnegative integers let 

si(w)^j2(-i) j ( d+1 ) (j + iyu^-i. 

3=0 

Lemma 11. The polynomial (w) has all positive real and simple roots. 
Proof. We have the following recursion (proof is straightforward) : 

si +1 (w) = (d + 2)s d i ( W )- W A- S f(w). 

Now apply [34( Vol. 2, Part 5, Problem 66] with necessary modifications in order to show that the successive 
application of the recurrence to Sq (w) = (w — 1) +1 does not introduce new multiple roots (but rather 
decreases the multiplicity of w = 1) . □ 

By applying Sturm's theorem on the interval (0, +oo) and using the fact that -^sf (w) = (d + 1) sf" 1 (w) 
we also get the following result. 

Corollary 12. All the polynomials sf (w) with i > d have positive real and simple roots. The sequence 
{ s d' s rf _1 > . . . ,s°} is a Sturm sequence. 

Proposition 13. Let u = zw (recall that z — u N ). Then the polynomial pfj (u) satisfies 

d 

(4.2) p d N (zw) = z d+2 Y^^N l s d (w). 

i=0 

Proof. Straightforward. □ 

Now let {ui = z, . . . , Ud] denote the roots of pfj (w), and {wi = 1, . . . , Wd} denote the roots of (w). 
Corollary 14. The distance between the {u\, . . . ,Ud} remains 0(1) as N — > oo. 

Proof. Consider the decomposition (|4.2[) . By Lemma [TTl {w±, . . . ,Wd} are positive, real and simple roots of 
(w). By Rouche's theorem, as N — > oo the roots of -j^zp% (zw) converge to {w\, . . . , Wd}- Obviously the 
polynomials jjzp% (zw) and p% (zw) have the same roots. □ 

Proposition 15. The derivative of pf^ (u) at U\ — z is equal to (d + l)\N d ctdZ d+1 . 

Proof. This is a combinatorial exercise. □ 



Proposition 16. The asymptotic estimate \-^P% (u) |«=u 4 1 ~ \ a d\ N d holds for all {ui}. 

Proof. Use the decomposition (|4.2[) . □ 

Now we can estimate the deviation of the roots of qfj from {ui, . . . , Ud}- 
Lemma 17. Denote by {yi, . . . ,yd} the roots of qfj. Then there exists Cg such that for N 3> 1 and for 

i = i,2,...,d 

\y j -u j \<C^N- d -\ 
Proof. By Proposition [111] there exists a constant Cio such that 

> C w BN d . 



-^p% (u) \ u=Ul 



Now again from (|4.2[) it is easy to see that the second derivative of p% at {ui, . . . , Ud} is bounded: for some 
constant Cn we have 



d 2 



< CnAA d . 



Take disks of radius rj(N) — Cg^N d 1 around each root U{ (where Cg is to be determined). By a 
straightforward application of Taylor's theorem (see e.g. [9l Lemma A. 6]) we therefore have that 



-l 



\p% {m + rj{N))\ > C 12 B V (N) N d = C 12 C 9 RN 

Now consider the perturbation polynomial e% A = qfj—pfj. Its coefficients have magnitudes \mjN — n^jN | < 
RN- 1 . Therefore 

\e% {ui + r,{N))\ < daRN- 1 . 

Note that for N 1 the constant C13 does not depend on Cg because, say, \ui + 17 (A)| < 2 \ui\ < (d). 

Consequently, if we choose Cg > ^| we can apply Rouche's theorem and conclude that qfj has a simple 
zero within distance Cg-giV _<i_1 from Uj. 

By Corollary [M] the {u^ are O (l)-separated, therefore if N is large enough than the quantity Cg^N^ d ^ 1 
will be smaller than the minimal separation distance. □ 

Remark 18. This analysis is valid for any root of qfj, not just the perturbation of u\ = z. The roots of qfj 
all lie approximately on the ray with angle £N. This means that the parameter £ can be recovered with high 
accuracy from any root of qjq (it), and we expect that it might be important for practice (so for instance one 
can approximate £ by averaging). 

Proof of Theorem^ first part. By Lemma H7] the accuracy of step 2 is bounded by Cg-giV~ d ~ 1 . Extraction 
of iV-th root in step 3 further increases the accuracy by the factor jj. As mentioned above, this operation 
introduces multiple-valuedness which can be circumvented by assuming that we know a-priori sufficiently 
close approximation to the true root £ (sufficiently close here should mean at least O (A -2 )). Step 4 preserves 
this estimate (by Taylor expanding of log (1 + e)). □ 



Proof of Theorem^ second part. We have recovered the approximate value ujn which satisfies |ojjv — < 
C 4 § N- d - 2 , while \fh k - m k \ < Rk" 1 . 
The error vector satisfies 

ag - a 
5i — a\ 

ad - otd 

It is immediate that 

V# = V 1 d diag{l,7V,...,iV d }. 



rriNU N — ijinlo 

~ 2N — IN 

m 2 N^> N — m 2N uj 



(d+l)N 



-(d+l)N 



Since 

\m jN \ < C u AN d 

= m jN + R(N)N- 1 , R(N)<R 

r ~ d - 2 , C 4 (N)<C 4 

then we have (using standard Taylor majorization techniques, see e.g. (9[ Proposition A. 7]) that 



R 

uj n = uj + d (N) —N~ 

B 



J N 



R 



-d-1 



B 



-jN 



Buj 



-jN 



R 



1 - Cib (N) -N 



-d-1 



and consequently 
(4.3) 



~ jN -jN 

rrijN^M — vrijjssui J 



■ C\,;R (l + ^l -V ' 



IjffW 


-jN 


and also Cn = f 


(vfy 1 


| . Since 




a - 


a 




"1 






"Co" 


Q.X - 


«i 




N- 1 






Ci 














a-d - 


Oid_ 













and using the estimate f|4.3[) we have immediately that for j = 0, 1, . . . , d 

A' 



\a j -a j \<C 5 R[l + -) N-*- 1 



where C5 = f CxqCh. This completes the proof of Theorem [51 □ 

5. Numerical experiments 

In our numerical experiments we compared the performance of the following Eckhoff-based methods for 
recovery of a single jump point position from the first M Fourier coefficients: original Eckhoff 's formulation 
from [17] (Eckhoff); our previous method from [9J (Half); the method presented in this paper (Full). All 
the three methods in essense solve a polynomial equation pm (u) — satisfied by the jump point u> = e~^: 
for Eckhoff and Half this polynomial is constructed from consecutive samples k = M — d — 1, . . . ,M 
(Algorithm [I} . while Full uses the decimated sequence k = N, 2N, . . . , M. The only difference between 
Eckhoff and Half is the degree of pm (u): the former uses the full smoothness d while the latter uses 

All calculations were done using Mathematica software with high-precision setting. The results, presented 
in Figure l5Tl on page[9l agree well with the theory - Full presents an improvement of ~ M~ d / 2 compared 
to Half, and improvement of order M~ d compared to Eckhoff. 



6. Discussion 

Eckhoff 's original conjecture related exclusively to his own method (for k = M — (d + 2) K, . . . , M) while 
assuming maximal possible reconstruction order. As stated by him in |17| it is false because in the general 
case there is no cancellation (as is evident also from our experiments). Eckhoff himself clearly knew about 







(a) Order=3 



(b) Order=4 



loan 





(c) Order=5 



Full 

(d) Order=7 




Figure 5.1. Fixed order reconstruction. Full represents the algorithm of this paper, BY 
2011 referes to the method of [9] while Eckhoff denotes the original method of Eckhoff 
from [T7]. 

the cancellation phenomena for d = 1 (as is evident from the explicit calculations in [17]), but he apparently 
failed to generalize it correctly for higher values of d. 

While both [9] and the present paper describe a real phenomena of accleration of convergence, we are 
yet to find a fundamental mathematical reason which could somehow explain these phenomena. We believe 
that such an explanation could be based on the description of the so-called "Prony map" - in particular the 
geometry of its image. We have started this investigation in [8), and we hope to get additional insights in 
that direction following the results of the present paper. In particular, we would like to determine whether 
applying efficient algorithms for solving Prony-type systems can eliminate the need for localization via 
convolution (as we noted in [5] the localization step is currently not easy to implement in practice). 

In a related line of investigation, it would be natural to attempt apply the decimation technique to other 
reconstruction problems in algebraic sampling, such as reconstruction from moments, finite rate of innovation 
and others (0 El EQl EH E3 [35]). 

There appears to be a natural way to generalize the ID algorithm to higher dimensions, by a technique 
similar to [6J. We plan to present results in this direction separately. 



2 These unfortunate shortcomings may be partially attributed to the fact that in his experiments he chose to reconstruct 
functions which were essentially either exact piecewise-polynomials or real analytic functions - cases in which the cancellation 
phenomenon always appears. 
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An "Eckhoff conjecture for the analytic case" (recall (|1.2[0 appears to be still open (but see a recent work 
[1]). One natural line of attack would be to analyze how the constants appearing in the accuracy estimates 
depend on the smoothness order d. 



Appendix A. Maximal accuracy for jumps 
Proof of Proposition^ Consider the following subset of the space PC (d, K) 



B{A,R)= |/G PC (d,K) : f = + \c k (*)| < R ■ k~ d ' 2 ; J2\a u \ < A 

where the smooth part *S> is in C d and the quantities {ai.j} denote the associated jump magnitudes of the 
piecewise polynomial whose degree is d. It is sufficient to show that there exists an absolute constant C 
such that for every index M there exist two functions gu , Hm £ B (A, R) whose first M Fourier coefficients 
coincide, while the corresponding jump locations differ by CM~ d ~ 2 . Once we show this, it is clear that 
no deterministic algorithm will be able to reconstruct the jump locations of all functions in B {A, R) with 
accuracy essentially better than O [M~ d ~ 2 ). 

Take gu — 9 to be a fixed piecewise polynomial with jumps {£1, ■ ■ • , Of} and associated jump 
magnitudes {aij}, satisfying Y^i j \ a i,j\ < A. Obviously g e B(A,R). Denote S — CM~ d ~ 2 where C is to 

be determined. Let $^ denote another piecewise polynomial of degree d with jumps 

{m = £1 + 5, . . . , r) K = £ K + 6} 
and the same jump magnitudes {ai j} as those of Let 

Ofe = Cfe 

Finally take 

fe=0 

Clearly c k (<?m) = c k (/im) for k = 0, 1, . . . , M. In order to ensure that Hm G B (A, R) we need to choose C 
small enough such that 

\a k \ < R- k~ d - 2 ; fc = l,2,...,M. 

def R 



$(<0 _ $w 

^ M 



M 



It turns out that C = # will fit just fine. Indeed: 



a k = c k 



if 



i+i 



3=1 



Now obviously 



if d 



a/. 



z=o ( ifc ) 



< 



2~7rfc' 



From geometric considerations we have 1 — e < 5 A;, therefore 

|« fc | < — 5k =—CM- d - 2 = — ■?-■ N-r d - 2 <RM- d - 2 . 
1 1 ~ 27rfc 2tt 2tt A 



This completes the proof. 



□ 
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